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1. INTRODUCTION 

In the last decades, stroke illness has been considered as one of the most well-known diseases for 
humans in the all sides of the world. Stroke illness has been considered as the main cause of human nervous 
system damage, lower limb function disorders and hemiplegia [1]. As aresult, the rehabilitation process of 
stroke patients has become an urgent thing to be undertaken to solve the disorder motion of aging humans that 
suffer from the stroke case [2]. Traditionally, the classical physical therapy is performed by rehabilitation 
therapists in manual manner. However, this manual process become tedious and exhausted strategy to the 
physical therapists while they try to help the stroke patients in process of recovering the gait of them. As a an 
effective and sufficient kinds of the rehabilitation robots, lower limb rehabilitation exoskeletons (LLE) have 
been used as an effective rehabilitation approach for stroke patients with motion disorder [3]. Various types of 
of LLE have been manufactured and constructed by different universities and research centers such as Lokomat 
[4], berkeley lower extremity exoskeleton (BLEEX) [5] and active leg exoskeleton (ALEX) [6]. Different 
control strategies have been introduced in order to track predeterminded trajectory for nonlinear robotics 
systems such as: active disturbance rejection control [7], fractional order control [8], robust control [9] and 
sliding mode control (SMC) [10]. 

In fact, SMC has been adopted as an effective and adequate controller for trajectory tracking of linear 
and highly coupled non-liner robotic rehabilitation systems because of having the distinqushed properties of 
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SMC represented by the insensitivity against the variations in parameters and the external perturbations [11], 
[12]. These unique aattributes enable the SMC to be used for diifrent applications with various combination 
such as: impedence control for rehabitation exoskeletons [13], and hybrid adaptive robust control for lower 
limb exoskeleton [14]. The sliding mode control is also used as an efficient controller for delta parallel robot 
[15] and for wind turbines [16]. The tuning parameters of SMC controllers represent an important issue on the 
system performance. The incoorect tuning of the parameters may lead to decreasing the controller performance 
[17]. So, to get an exccellent performance, a meta-heuristic optimization technique adopted from natural 
behavior of animals has been used in last decades. In the present study, grey wolf optimization (GWO) has 
been used to optimize the tuning factors of developed sliding mode controller. Traditionally, the dynamic 
equations of lower limb robotic exoskeletons have been developed using popular Lagrange equations [18]. In 
fact, the contacting of robotic systems with horizontal or vertical surface modified the Lagrange equations to 
include the Lagrange multiplier and Jacobian vectors in its formula [19], [20]. Most of previous works in 
literature have modeled LLE without taking in to consideration the contacting of LLE with the ground 
(constraint motion). Hence, the primary objective of this article is to model and control the LLE in case of 
contacting with ground. To the best of author's knowledge, the modeling and controlling of LLE in constraint 
motions have not been discussed deeply in literature. The rest of this paper is constructed as follows: The 
derivation of dynamic equations of the LLE in constrains motions is accomplished in the section 2. The optimal 
sliding mode control law are developed in section 3. In the section 4, the results of the work have been explained 
and discussed. The conclusion of this article is achived in the section 5. 


2. METHOD 
2.1. Modeling of lower limb exoskeleton in free motion 

The dynamic model of the LLE with three degrees of freedom has been adopted in this study as indicated 
in Figure 1. Figure 1(a) represents the free motion of LLE, whereas Figure 1(b) stands for constrained motion of 
LLE. In free motion, classical Lagrange equations are used to model the dynamic equation of LLE as shown 
Figure 1(a) where O is the coordinate origin and h refers to the distance from the coordinate origin O at hip joint 
to the point of contacting of LLE with the ground and @(q) refers to an algebraic constraint equation that represent 
the constraint motion in joint space [21], [22]. The Euler-Lagrange formulas are normally used to derive the 
dynamic model equations of the lower limb exoskeleton, which can be expressed as: 


M(q)g+C(q,q)qt+G(q)+tq=T (1) 


$(9,.92-9;) #9 


, Y 
y Ground Ground 
(© Hip joint © Knee joint © Ankle joint © Hip joint Ce) Knee joint © ankle joint 
(a) (b) 


Figure 1. Three-link exoskeleton structure in (a) free motion and (b) constraint motion 


where: q,q,q@ € R?are angular joint position, velocity and acceleration vectors respectively. M(q) € 
R3*3 stands for positive definite inertia matrix. C(q, q)q € R? represnts Carioles, centrifugal forces and torques. 
G(q) € R® are the torques of the gravity and t € R? is the vector of joint torques and tg € R? is the bounded 
unknown external disturbances at LLE joints. The contents of M(q) , C(q, q)q and G(q)of three joints of LLE 
are: 


M,, =a, + a, +a, + 2.a3cos(qz) + 2a. cos(qz + G3) + 2a,cos (q3) 
My2 = M2, = az + a4 + A3C0S(q2) + As cos(qz + gz) + 2a¢cos (q3) 
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M,3 = M3, = a, + as cos(qz + G3) + 26 cos(q3), Moz = az + Ay + 2a, cos(qs) 
M23 = M304 + 2a, cos(q3), M33 = ay 
Cy, = —As(q2 + G3) sin(q2 + G3) — A3q2 sin(q2) — agq3sin (q3) (2) 
Cyp = Os (ds + Go + ds) sin(4z + 3) — a3 (4s + G2) sin(qz) — agdgsin (aa) 
Cy3 = —a5(qi + dz + G3) sin(qz + G3) — a3 sin(qz2) — a6(qi + do + q3)sin (q3) 
Cy1 = —as5q,sin(qz + q3) + a3q; sin(qz) — a¢q3 sin(q3) ,C22 = —agq3sin (q3) 
C23 = —A6(Gi + G2 + Gz)sin (q3) 


C31 = Aq, Sin(qz + G3) + A6(G1 + G2) sin(q3), C32 = G6 (Gi + G2) Sin(q3) , C33 = 0 


G, = —[b,sin(q,) + bzsin(q, + gz) + bssin(q; + q2 + 43)] 
Gz = —[bzsin(q, + q2) + bssin(q, + q2 + q3)] ,G3 = —[b3sin(q; + q2 + 43)] 
by = (mM, 44 4Mz1,,M3l,)g , bz = (Mzd24Mzl2)g , bz; = M3d3g 


ay = iL, + m,d? + (m, + m3)l? Py az = IL, + mdi + m3 5 a3 = (m2d, + M3l,)l, 
a, =1,+m3d2 as = mgdgl, ag = mad3l, 


where: m;m,zm31,1,1, refer to masses and lengths of thigh, shank, and foot links of the human lower limb 
and exoskeleton respectively. d, dz dz stand for the position of the canter of mass of thigh, shank and foot 
of LLE and human lower limb respectively. 1,121, Stands for the moments of inertia of thigh, shank and 
foot of the exoskeleton and human lower limb respectively and g is the gravity acceleration. 


2.2. Modeling of lower limb exoskeleton in constrained motion 

During constrained motion, the LLE is in contact with ground and can be viewed as a closed loop 
chain. Thus, as shown in (1) of free motion can not be applied when we want to describe the dynamic 
equations of the LLE when contacting with the ground as shown in Figure 1(b). In this case, a holonomic 
constraint illustrated by an algebraic equation in joint space is used to describe the contacting of LLE with 
ground [22], [23]. 


0(q) = l, cos(q,) + l, cos(q, + qz) + 13 cos(qi + dz + q3) —h (3) 


Owing to theses constraint circumstances, the term t = J(q)"A must be added to dynamic in (1). The Lagrange 
multiplier A refer to the forces of contact of the LLE with the ground. Hence, the constrained dynamic equation 
of LLE can be written as (4) [19]: 


M(q)G+C@,4)4 + G(q) + tq =t+J(q)'A (4) 
where J(q) = ae stands for the Jacobian of the equation of constraint with: 
Ca ee ap is 
= Oe sin(q,) — lz sin(q; + q2) — ls sin(qi + q2 + 93) 
49(q) , ‘ ag : 
i> oo = —l, sin(q, + q2) — ls sin(q, + q2 + 43) J3 = oe = —l; sin(q, + q2 + 43) 


This constraint motion makes LLE lose number of degrees of freedoms which are equal to number of 
constraints [19]. In our work, the LLE is contacting with ground and this thing will make it lose one degree of 
freedom because of participating of one Algebraic equation. This thing separates the number of the degree of 
freedom of LLE into two main groups: the first group includes the independent joints coordinates q, = 
[q1,492]", 4, € R?.whereas the second group includes dependent joint coordinates. gy € R+, qy = qz . Hence, 
, areduced order model has been obtained that relay on independent joints only (Hip and Knee joints) and the 
motion of ankle joint will depend on the independents joint motion. If we differentiate the algebraic constraint 
in (3) with respect to time, we obtain: 
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p d@ 00(q) d ; A , : 
0(q) = ae rare =Jq =idi + J2d2 + J3d3] = 0 (5) 


the dependent velocity q3 can be expreesed in terms of independent velocities q, , dz as (6). 


ds = —J3(q)*h@) —Js(Q)“J2(@) (6) 
Hence, velocity coordinates ( 41, q2, q3) can be expressed as (7): 
[ o | ; i | | (7) 
dsl I“) *h@) = —Js(a) 2a) 
thus, we can get (8): 
q = H(q)4, (8) 


whereg=[41 42 4]".g,=[% G2)’, H = 


1 O ; / 
_ J1(q) _ la@ 

0 1 , Jaa = Jz oo 
J3(q) J3(q) 


Z1 Jz2 


The generalized acceleration can be written as (9): 


@ = H(q)q, + H(q)4z (9) 


if we Sub (9) in (4), the dynamic in (4) can be expressed as (10): 


M(Q)H(Q) Gz + C4, 9)d2 + GQ) + ta =T+J(Q)'A (10) 
where 
C,4,.4) = C@,qH@) +M@)H@) (11) 
with 
pe ; : | de = _ JaqJu@li@ls@ ins Ja(qJa(qy—J2@Jac@) 
i 172% Ua(q))2 122 Ua(q))2 
and 


Jiqa=¥1h + Wod2 + W3q3, Jaq=¥2h + W2d2 + W343, JI3q)=¥3 + W3d2 + W343 
W, = —[L, cos(q1) + lz cos(qy + G2) + l3 cos(qi + qz2 + 93)] 

W2 = —[l2 cos(qy + qz) + lz cos(qi + 92 + 93)] 

3 = —[ ls cos(qi + G2 + G3)] 


According to (10) refrs to the reduced order dynamic model for LLE when it contacting with ground. This 
equation includes J(q)1A term in its structure. The deleting of this term is very essential for developing an 
efficient controller. The J(q)™A term from (10) can be eliminated, by using the expression J(q)™H"(q) = 0 
which is popular property in literature which is used for different fields [22]. If we borrow this property and 
pre multiply equation (10) by H7, we can get (12) [23]: 


M(q)Gz + C(q,q)dz + G(q) + Tq =T (12) 
With M(q) = H™M(q)H ,C(q,q) = H"C, (4,4), G(q) = H"G(q), t = H"1, Tq = H"Tq 


3. CONTROLLERS DEVELOPMENTS 

In what follows, a sliding mode control scheme (SMC) will be developed in subsection 3.1 which 
assumes that the tuning of parmeters is peformed manually by trial and error strategy. Then, in subsection 3.2 
a grey wolf optimizer (GWO) will be used in order to get an optimized version of the developed controller 
(GW-SMC). Later, the description of grey wolf optimizer will be explained in detail in subsection 3.3. 
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3.1. Sliding mode control (SMC) 
Firstly, we present the following sliding manifold: 


S, = é, + Ae, (13) 


with 22-42 ~ qe ed, —G2 S, = [S21,5,2]"and A € R?*? , are positive diagonal matrices. The subscript z 


stands for independent joint coordinates which are the Hip and Knee joints in our paper. By taking the 
derivative of (13), one obtains: 


S, = 6, + Aé, (14) 


it is essential to enforce the system to be in the sliding surface. Hence, we obtain: S, = 0, then we obtain: 


de — 92 + ACG. — a) = 0 (15) 
by re-writing (12) and dropping dependency, we can get: 

G, = M[t— Cq, — G] (16) 
sub (16) in (15), we can obtain: 
M*[z — Cq, — G] — a2 + Ad, — af) = 0 (17) 
if we multiply equation by™ , we get: 

[= — Cq, — G] — Maf + MAG, — a8) (18) 


consequently, we can obtain the following sliding mode control law: 


Pat Pie (19) 
with 

Tequ = Cqz + G + Mg? — MA(G, — 42) (20) 

Tsy = —MK,sgn(S,) (21) 


where: T,,, is a high frequency discontinuous term which has been added to guarantee robustness with K, € 
R?*? is diagonal positive matrix. Thus, we can get the following sliding mode control law. 


Tsuc = Cdz + G + Ma? — MA(q, — a¢) —MK,sgn(Sz) (22) 
If we sub (22) in (14) we get: 

S.=f.,M* (23) 
the lyapunov function related to dynamic system of LLE is: 


v=i 


~ 2 


T 
Sz Sz (24) 
the main aim of this selection is ensure the minimization process of S, = 0. In addition, to remain on sliding 
surface and to ensure that the the error e, will convergence to zero. Hence, if we differntate V with respect to 
time, we get: 
Vass, (25) 
Vea S53 tag (26) 
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V = —S,"K,sgn(S,) < 0 (27) 


One of the main drawbacks of the sliding mode control is the existing of chattering problem which may cause 
the damage of the actuators. For eliminating and overcoming the influence of chattering effects, we replace the 
sign function by a smooth hyperbolic tangent function which has [24]. 


eSz—e-Sz 


tanh(S,) = (28) 


eSz4+e—Sz 

3.2. Optimal sliding mode control (GW-SMC) 
As shown in (22), the tuning of controller parameters of SMC ( A, K, ) has been usually performed 

using trial and error method. In fact, this is an exhausted process. For solving this tediedous method, a meta- 

heuristic optimization technique called grey wolf optimizer (GWO) presented in detail in sub-section 3.3 is 

used to tune the gains of SMC controller. Hence, we get aGW-SMC controller which can be expressed as: 


Tow-smc = Cd, + G + Ma? — MA*(q, — q¢) —M,K,*sgn(Sz) (29) 


where: ( A*, K,"), represent optimized values of SMC parameters. The block diagram of SMC optimized by 
grey wolf optimizer (GWO) is shown in Figure 2. 


3.3. Grey wolf optimizer algorithm (GWO) 

A grey wolf optimizer algorithm (GWO) is arelatively new swarm population-based optimization 
algorithm developed in 2014 by Faris et al. [25]. The work of GWO mainly depends on simulateing the social 
behavior of grey wolves that lived in the wild by considering the concept of the eadership hierarchy and the 
mechanism hunting to the prey of grey wolf. In this algorithm, a grey wolf hierarchy has been classified into 
four types of wolves: Alpha (a) wolf, Beta (8) wolf, Delta (6) wolf and Omega (@) wolf. The hunting process 
decision is achived by leaders of whole group which is Alpha (a) wolf. Beta (B) wolf follows the commands 
of Alpha (a) wolf and help it in decision of hunting and other activities in the pack. Omega (@) wolf will lead 
the other Omegas (@) and will follow Alphas (a), and Betas (8) wolves. If the wolf it is not an Alpha (a), Beta 
(B), or Omega (), the wolf will named as Delta (5) wolf. The optimal solution which is representing the 
location of prey will be implemented by estimating the populations of wolves by using an iterative method. 
The best wolf is an Alpha (a) wolf, the best second solution stands for Beta (8) and Delta (6) is third-best 
solution whereas Omega () wolves are the least significant solutions. Equations (30) and (31) have been used 
for formulating a mathematical model that describes the behavior of wolves in encircling process (finding the 
optimum solutions): 


D =|Cx,(t) - X()| (30) 
X(t +1) = |X,(¢) — AD| (31) 


where t refer to current iteration, A and C are coefficient vectors, Xs (t) stands for the prey victim position 


vectors. Whereas X (t)is the position vector of a grey wolf. The calculation of vectors Aand C can be expressed 
as: 


A=2ar%,-a (32) 
C=27 (33) 


1, Tz stands for vectors in random with range between 0 and 1, a decreased from 2 to 0 in linear manner during 
the iterations process. The updating process can be expressed as: 


Da = |C,Xq(t) — X|, Dg = |C.Xg(t) — X|, Ds = |CzX5(t) — X| (34) 
£,-= |%p= A,Da| X= [Ky =ADs| /Xq = [Xp —AsDs| (35) 
¥e41) = tess (36) 


3 
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the details of Pseudo code of GWO can be found in [25]. An integral time absolute error (ITAE) expressed in 


(37) has been used as an objective function to evaluate the position of each serch agents (wolfs) during the 
search for choosing the best value [17]. 


ITAE = f, tle| dt (37) 


fF] 


Figure 2. Block digram of sliding mode controller optimized by GWO (GW-SMC) 


4. RESULTS AND DISSCUSION 

This section is dedicated for comparing the effectiveness of developed optimal sliding mode controller 
GW-SMC derived in the previous section with recentely developed SMC with chateering suppressor which is 
develpbed in [26]. The desired trajectories are taken from [27]. The parameters of human and exoskeleton that 
used in simulation are adopted from [28] and listed in Table 1. Moreover, the integral absolute value error 


(IAE) and integrated squared error ISE have been used as a measure for tracking performances of SMC and 
GW-SMC controllers. 


IAE = fle| dt (38) 

ISE= fete dt (39) 

The parameters values SMC controller with with chateering suppressor has been selected by trial and 
error method. For GW-SMC controllers, the parameters selected as A = 25.73, K, = 8.22 and have been 
tuned using GWO. The number of wolves chosen is 50 and the algorithms are repeated for 70 iterations. Figure 


3 explain the excellent Objective Function performance of GWO after 70 iterations. The cpmparative 
simulation results are explained for three cases: nominal, uncertainty and disturbance rejection tests. 


Table 1. Physical parameters of human and lower limb exoskeleton [28] 


parameter Mass (m) in kg Length (1) in m Inertia (Thin kg.m2 Center of length (d) m 
Exoskeleton Human Exoskeleton =Human _ Exoskeleton Human Exoskeleton — Human 

Thigh 0.2043 7.33 0.41 0.407 5.7x10-3 0.1502 0.15 0.1763 
Shank 0.2159 3.4503 0.39 0.4334 4.3x10-3 0.0505 0.11 0.1849 
Foot 0.115 1.075 0.159 0.275 3.9x10-4 0.0038 0.04 0.1179 


4.1. Nominal case 

The evolutions of position, position error, and torque input signal for hip and knee joints for SMC in 
[26] and GW-SMC controllers in nominal case are depicted in Figures 4 and 5 respectivaly. It is worth to 
observe the validity of the GW-SMC controller over SMC controller with chateering suppressor in [26] in 
terms of trajectory tracking performance and minimization the tracking errors. This refelects the vital role of 
GWO in finding the optimal values and its impact to the system performance. A smooth control input signals 
have been obtained for two controllers owing to existing hyperbolic function which has the ability to reduce 
the chattering effect effectively. The results of IAE index and ISE performance indices for nominal cases for 
SMC in [26] and GW-SMC controllers are shown in Table 2. The results explain that GW-SMC controller has 
lowst tracking error as compared with SMC controller with chateering suppressor. 
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Figure 3. Evolution of objective function of grey wolf optimizer after 70 iterations 
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Figure 4. Performance of hip joint for two controllers in nominal case 
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Figure 5. Performance of Knee joint for two controllers in nominal case 
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Table 2. Performance indices IAE and ISE for two controllers for nominal case 


Index 


IAE ISE 
Controller GW-SMC (29) SMC [26] GW-SMC (29) SMC [26] 
Hip joint 0.001372 0.013715 0.001654 0.01571 
Knee joint 0.001541 0.005277 0.006407 0.01068 
Both joints 0.002914 0.018992 0.008062 0.02640 


4.2. Uncertainty case 
In this case, the parameters of the LLE 


are varied by 30% from their original nominal values. It is 


observed from Figures 6 and 7 and Table 3 that GW-SMC controller has an excellent robustness against system 
uncertainty. The figures and table illustrate that the performance of SMC with with chateering suppressor is 
significantly changed in case of the presence of parameters variations while the proposed GW-SMC is remain 


insensitive to parameter variations. 


Position of hip joint (rad) 


desired 


GW-SMC SMC[26] 


— — smci26]| 
rh 


ee ae aieiae ee 


time(sec) 


Figure 6. Performance of hip joint for two controllers in uncertainty case 


Position of Knee joint (rad) 


desired GW-SMC SMC[26] 


wetcctccee =] 


time(sec) 


Figure 7. Performance of Knee joint for two controllers in uncertainty case 


Table 3. Performance indices IAE and ISE for two controllers for uncertainty case 


Index 


IAE ISE 
Controller GW-SMC (29) SMC [26] GW-SMC (29) SMC [26] 
Hip joint 0.001649 0.021897 0.001933 0.02523 
Knee joint 0.001705 0.013445 0.006690 0.01962 
Both joints 0.003355 0.035342 0.008623 0.04486 
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4.3. Disturbance rejection case 
The third simulation has been achieved in the presence of time varying disturbance which has the 
following expression [29]. 
: a ? sin(t) + 0.5sin (100¢) (40) 
4 =|Tg2) =| cos(t) + 0.5 sin(100t) 


We can notice from Figures 8 and 9 that smaller control errors and faster convergence have been effectively 
ensured by GW-SMC as compared with the SMC with chateering suppressor in [26]. Thanks to tuning process 
by GWO our GW-SMC controller can still ensure the best comprehensive control performance and it is still 
satisfactory in spite of the presence external disturbances. The tracking control performance indices [AE and 
ISE for this case are given in Table 4. We can notice that the GW-SMC has the lower best values in terms of 
IAE and ISE performance indices as compared with SMC with chateering suppressor. To sum up, the 
advantages and the superiority performance of the proposed GW-SMC controller have been verified. GW- 
SMC controller provides an excellent tracking performance; high precise tracking, less chattering and good 
robustness towards the uncertainties and external disturbances. 


Position of hip joint (rad) 
T 


T desired GW-SMC SMC[26] 


time(sec) 


Figure 8. Performance of hip joint for two controllers in disturbance rejection case 


Position of Knee joint (rad) 


desired GW-SMC SMC[26] 


Position Tracking error of Knee joint (rad) 


time(sec) 


Figure 9. Performance of Knee joint for two controllers in disturbance rejection case 
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Table 4. Performance indices IAE and ISE for two controllers for disturbance rejection case 


Index JAE ISE 
Controller GW-SMC (29) SMC [26] GW-SMC (29) SMC [26] 
Hip joint 0.0026758 0.026384 0.003066 0.032597 
Knee joint 0.002160 0.046239 0.007396 0.048754 
Both joints 0.004836 0.072624 0.010463 0.081352 


5. CONCLUSION 

In this study, a reduced order constrainnt model with augmented Lagrange equations for a lower limb 
exoskeleton in constrained circumstances dedicated for rehabilitation of stroke patients have been developed. 
Optimal sliding mode controller tuned by grey wolf optimization algorithm (GW-SMC) has been developed to 
control lower limb exoskeleton in constrained enviroinments. The stability analysis of the closed-loop system 
has been performed using Lyapunov theory of stability. Simulation results in nominal, uncertain and 
disturbance rejection tests indicate that the sliding mode controller optimized by grey wolf optimizer provides 
excellent results performances in terms of robustness, fast response and tracking error as compared with 
conventional sliding mode controller. In the Future, the focusing will be on extending this controller in real 
time rehabilitation exoskeleton systems and using another types of optimization algorithm for tuing the 
controller parmeters. We are planning to combine this controller with fractional ordr calculas for getting an 
excellent response of the performance of the system. 
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